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IMTRODUCTIOH 


The  development  of  computational  fluid  dynamics  (CFD) 
procedures  has  progressed  extremely  rapidly  during  the  past  two 
decades.  The  parallel  rapid  development  in  computer  hardware 
resources  and  architectures  has  not  only  matched  the  explosive 
algorithm  development  but  has  indeed  provided  and  continues  to 
provide  its  impetus.  Together,  the  resources  are  now  available  for 
the  numerical  simulation  of  the  flow  about  complex  three 
dimensional  aerospace  configurations. 

A  particularly  large  growth  has  occurred  during  the  last  few 
years  in  the  computer  hardware  sector  of  scientific  workstations. 
Many  vendors  are  supplying  fairly  advanced  systems  for  under  twenty 
thousand  dollars.  Is  such  a  workstation  sufficient  for  solving  for 
the  compressible  viscous  flow  about  a  complete  aerospace 
configuration,  perhaps  representing  a  leveling  of  the  playing  field 
previously  reserved  for  those  with  access  to  supercomputers? 
Memory  capabilities  appear  to  be  expanding  rapidly,  faster 
apparently  than  comparable  improvements  in  hardware  speed.  The 
efficiency  of  the  new  numerical  algorithms  may  be  the  deciding 
factor  in  answering  the  above  question,  at  least  for  the  near 
future . 

In  reality,  the  workstation  will  not  replace  the 
supercomputer.  However,  many  of  the  tasks  now  performed  on 
supercomputers,  such  as  the  important  task  of  program  development, 
can  be  down  shifted  to  workstations.  The  tasks  of  large  number 
crunching  can  be  left  to  the  more  powerful  supercomputers.  In 
computational  fluid  dynamics,  program  development  requires  many 
runs  of  a  developing  computer  code.  In  three  dimensions,  even  on 
coarse  grids,  these  runs  can  be  long  enough  to  require  the  use  of 
a  supercomputer.  Efficient  algorithm  development  toward  faster 
flow  solvers  is  shortening  the  durations  of  these  test  runs  so  that 
in  time  more  and  more  of  them  can  be  made  on  inexpensive 
workstations . 

Hypersonic  flow  has  been  an  area  of  much  scientific  interest 
in  recent  years.  The  description  of  a  gas  flowing  about  a  body 
traveling  at  hypersonic  speeds  is  very  complex.  It  must  include 
finite  rate  chemical  reactions  as  the  freestream  gas  dissociates 
upon  passing  through  the  very  strong  bow  shock  wave.  Also,  because 
these  changes  of  state  occur  so  rapidly,  the  gas  can  be  expected  to 
be  in  thermal  nonequilibrium,  requiring  the  determination  of 
separate  translational,  rotational  and  vibrational  temperatures  of 
the  gas.  These  two  additional  gas  features,  chemical  reaction  and 
thermal  nonequilibrium,  require  the  solution  of  many  additional 
equations.  A  perfect  gas  description  in  three  dimensions  requires 
one  continuity  equation,  three  momentum  equations  and  a  single 
energy  equation.  On  the  other  hand,  a  nonequilibrium  gas  requires 
a  continuity  equation  for  each  species  present,  three  momentum 
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equations  plus  an  energy  equation  for  each  thermal  mode  present. 
For  a  seven  species  model  for  air  (N2,02,N0,N0*,N,0,e')  and  a  three 
temperature  gas  model  (T,Tr,Tvib) ,  thirteen  equations  are  required 
as  compared  to  five  for  a  perfect  gas.  Because  the  gas  description 
is  very  much  more  complicated  for  hypersonic  flow,  the  efficiency 
of  the  numerical  solution  method  is  of  high  importance. 

The  main  purpose  of  the  report  is  to  present  an  efficient 
numerical  method  for  solving  the  equations  of  viscous  compressible 
flow  in  three  dimensions.  A  main  goal  of  the  present  study  was  to 
apply  this  method  to  simulate  the  flow  about  a  generic  hypersonic 
vehicle  in  powered  flight.  This  is  a  rather  monumental  task  that 
requires  very  high  numerical  efficiency  to  make  the  calculation 
possible.  The  method  presented  herein  is  fully  implicit  and  uses 
block  tridiagonal  line  inversion  with  Gauss-Seidel  relaxation.  The 
block  tridiagonal  procedure  requires  the  use  of  a  structured  grid. 
This  type  of  algorithm  has  been  very  successful  in  two  dimensions 
where  the  line  direction  extended  from  the  body,  through  the 
boimdary  layer  fine  grid,  across  the  shock  wave  and  into  the 
freestream  and  where  the  streamwise  direction  was  treated  in  a 
Gauss-Seidel  manner.  This  numerical  procedure  typically  converges 
for  steady-state  solutions  in  under  100-time  steps  instead  of  the 
usual  several  hundred  or  thousands  for  viscous  compressible  flow 
simulations.  The  extension  to  three-dimensions  has  proved  to  be 
difficult.  To  avoid  a  directional  bias,  a  block  tridiagonal  line 
inversion  procedure  is  required  in  the  spanwise  direction  as  well. 
The  stitching  together  of  two-line  inversion  procedures  with  Gauss- 
Seidel  relaxation  to  form  a  numerical  method  of  comparable 
efficiency  to  the  two-dimensional  method  will  be  presented. 

A  computer  program  based  upon  this  method  has  been  coded.  The 
computer  code  has  not  been  vectorized.  It  has  been  run  only  on  a 
serial  processing  workstation.  In  addition,  for  developmental 
purposes,  the  program  has  been  witten  so  that  Jacobians,  for 
example,  have  been  written  as  subroutines  and  called  repeatedly  by 
other  parts  of  code.  This  is  less  efficient  than  placing  them 
where  needed  within  the  code,  but,  on  the  other  hand,  modifications 
to  the  code  need  to  be  made  only  at  one  place.  Furthermore,  the 
Jacobians  of  the  flux  vectors  in  each  coordinate  direction  (SF/SU, 
SG/SV  and  SH/SG  etc.)  have  been  inritten  as  a  single  generalized 
Jacobian  with  rotational  direction  cosines  that  can  be  used  to 
determine  each  individual  coordinate  direction  Jacobian.  Again, 
Jacobian  code  modifications  using  this  formulation  require  changes 
in  only  one  place  in  the  computer  program.  For  the  computational 
efficiency  of  production  runs,  unlike  code  developmental  runs,  this 
formulation  should  be  removed  and  the  resulting  code  should  be 
vectorized  to  the  fullest  degree  possible. 

The  computer  code  was  applied  to  two  test-flow  problems  for 
which  experimental  or  flight  data  was  available  for  program 
validation.  The  first  case  simulated  the  experiment  of  Lobb  [1], 
1964,  for  Mach  15.3  flow  past  a  sphere  in  air  which  exhibited 
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chemical  nonequilibrium  and  thermal  vibrational  nonequilibrium 
behavior.  The  second  test  case  simulated  the  RAM-C  1 1  fliqht 
measurements  [2]  made  in  the  late  I960's  of  peak  electron  number 
densities  in  the  Mach  25.9  flow  past  a  sphere  cone  body  at  71  km. 
altitude.  This  flow  exhibited  ionization  and  translational, 
rotational  and  vibrational  nonequilibrium.  Finally,  the  computer 
program  was  applied  to  simulate  the  flow  about  a  generic  hypersonic 
vehicle  at  Mach  25  at  high  altitude.  The  body  consisted  of  a 
sphere-cone-cylinder  fuselage  with  a  delta  wing  attached.  The 
delta  wing  tips  were  folded  up  to  resemble  a  hypersonic  glider.  An 
engine  was  then  placed  below  and  attached  to  the  fuselage-wing 
surface.  The  underside  of  the  body  downstream  of  the  engine  exit 
was  contoured  to  form  a  nozzle  surface.  The  flow  entering  the 
engine  inlet  plane  was  accelerated  to  produce  thrust  at  the  exit 
plane.  Body  drag,  engine  ram  drag,  engine  thrust  and  surface  heat 
transfer  were  calculated. 

The  present  study  combines  the  PhD  thesis  research  on 
hypersonic  nonequilibrium  flow  of  three  former  Stanford  University 
graduate  students.  Dr.  Graham  V.  Cemdler  (1988)  [3],  Dr.  Tahir 
Gokcen  (1989)  [4],  and  Dr.  Xiao-lin  Zhong  (1991)  [5]  with  the  three 
dimensional  algorithm  research  presented  in  Refs.  [6,7]  to  study 
three-dimensional  nonequilibrium  hypersonic  viscous  flow. 


GOVERNING  EQUATIONS 

The  Equations  governing  nonequilibrium  viscous  compressible 
flow  can  be  written  in  vector  form  as: 


dt  dx  dy  dz 


where  U=[pi,p2, . . .  ,p„,pu,pv,pw,ei, . . .  ,e„]^,  F,G,  and  H  are  flux 
vectors,  and  W  is  a  source  term  vector.  The  flux  vectors  contain 
terms  representing  the  convection  of  species  mass,  momentum  and 
energy,  pressure  and  viscous  stress,  work  and  heat  transfer.  The 
viscous  and  heat  transfer  coefficients  are  obtained  from  Blottner,s 
curve  fits  [8]  and  use  of  Eucken's  [9]  and  Wilke's  relations  [10]. 
The  source  vector  W  contains  terms  for  the  production  and  loss  of 
specie  densities  through  chemical  reaction  and  energy  transfer 
through  thermal  relaxation.  Vibrational  relaxation  is  modeled 
using  Park's  modified  Landau-Teller  rate  equation  [3]  and 
rotational  relaxation  is  modeled  using  Parker's  formulation  for 
relaxation  times  [11].  Expressions  relating  vibrational  energies 
to  vibrational  temperature  and  energies  contained  in  excited 
electron  states  were  taken  from  Candler  [3]. 

The  air  chemistry  model  as  obtained  from  Candler's  PhD  thesis 
[3]  is  given  below  (see  top  of  next  page). 
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NUMERICAL  METHOD 

Chakravarthy  in  1984  [12]  presented  unfactored  implicit 
relaxation  methods  for  solving  the  Euler  equations.  By  avoiding 
the  technique  of  approximate  factoring  and  the  error  it  introduces 
that  slows  convergence,  his  approach  is  very  efficient  for  solving 
the  equations  of  compressible  fluid  flow.  The  approach  uses  type 
dependent  differencing  and  line  relaxation.  These  are  essentially 
the  same  techniques  used  earlier  by  Murman  and  Cole  [13]  for 
solving  the  transonic  small  disturbance  equation  and  later  in 
rotated  form  by  Jameson  [14]  for  the  full  potential  equation.  The 
type  dependent  differencing  for  the  Euler  equations  can  use  "flux 
vector  splitting"  introduced  by  Steger  and  Warming  [15]  or  now  the 
more  popular  "flux  difference  vector  splitting"  of  Roe  [16].  The 
line  relaxation  consists  of  block  tridiagonal  matrix  inversions  in 
the  across  the  flow  direction  with  Gauss-Seidel  relaxation  used  in 
the  streamwise  direction.  The  type  dependent  differencing  is 
upwind  and  creates  implicit  block  matrix  equations  with  stronger 
weights  given  to  the  diagonal  matrix  elements  than  would  otherwise 
occur  using  central  difference  approximations.  This  facilitates 
the  use  of  an  unfactored  relaxation  algorithm  for  the  solution  of 
the  matrix  equation  Instead  of  the  more  common  direct  approximate 
factorization  algorithm  or  direct  Gaussian  elimination  with  their 
inherent  numerical  inefficiencies. 

In  1985  efficient  Gauss-Seidel  line  relaxation  algorithms  for 
solving  the  Navier-Stokes  equations  in  two  dimensions  were 
presented  [17,18].  Their  accuracy  was  improved  for  calculating 
leuninar  and  turbulent  shear  layers  in  1987  [19].  Accurate  two- 
dimensional  Navier-Stokes  calculations  were  obtained  in  100-time 
steps  or  less  with  CFL  numbers  reaching  one  billion.  Also  in  1987, 
Candler  extended  the  technique  to  three  dimensions  for  flow  past 
aircraft-like  configurations  [20].  Block  tridiagonal  line  inversion 
was  used  in  the  normal  to  the  body  direction  and  Gauss-Seidel 
relaxation  in  the  streamwise  and  around  the  body  directions.  In 
1988  difficulties  with  this  method  were  encountered  for  three- 
dimensional  cascade  flow  [6].  Flows  that  should  remain  two- 
dimensional  became  three-dimensional  because  the  Gauss-Seidel 
procedure  of  using  the  latest  available  data  introduced  an 
asymmetry  as  each  plane  was  updated  using  new  data  on  one  side  and 
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older  data  on  the  other  side  of  the  block  tridiaqonal  "line".  In 
addition,  the  cascade  flow  required  block  tridiaqonal  inversion 
procedures  in  both  cross  stream  directions,  blade  to  blade  and  side 
wall  to  side  wall.  The  construction  of  a  three-dimensional 
unfactored  algorithm  containing  block  tridiagonal  inversion  in  two 
directions  is  not  trivial.  The  statement  found  in  many  algorithm 
development  papers,  including  this  author's,  "the  extension  to 
three  dimensions  is  straight  forward"  is  untrue. 

The  three-dimensional  implicit  algorithm  used  for  the  results 
presented  later  does  combine  block  tridiagonal  inversion  in  two 
directions  with  Gauss-Seidel  relaxation  in  the  third  direction.  It 
is  first  order  accurate  in  time  and  second  or  third  order  accurate 
in  space.  It  could  be  made  second  order  in  time  also,  but  for  the 
problems  considered  here,  for  which  rapid  convergence  to  steady 
state  was  the  goal,  order  of  accuracy  in  time  loses  its  meaning. 
The  governing  Navier-Stokes  equations  with  nonequilibrium  effects 
were  approximated  using  low  dissipation  flux  split  approximations 
for  the  Euler  terms  as  described  in  Ref.  [19]  and  central 
difference  approximations  for  the  viscous  terms. 

Consider  the  flux  at  surface  i+l/2  lying  midway  between  mesh 
points  (i,j,k)  and  (i+l,j,k).  The  flux  F  can  be  split  as  follows. 

where  A  is  the  Jacobian  5F/4U,  A*  and  A'  represent  the  split 
Jacobians  containing  the  positive  and  negative  eigenvalues, 
respectively,  of  A,  and  A  «  A*  +  A".  Ut.  represents  the  value  of  U 
at  the  surface  as  approached  from  the  left,  and  Ur  is  similarly 
defined  using  data  from  the  right.  At  flux  surfaces  that  are 
supersonic  on  both  sides  the  flux  is  not  split  at  all;  either  A*  or 
A"  is  null.  At  surfaces  within  subsonic  regions  of  the  flow  and 
away  from  shock  waves  third  order  flux  splitting  is  achieved  by  (1) 
using  an  interpolation  formula  and  the  nearest  four  points,  (i- 
lrj»h),  (i,j,k),  (i+l,j,k)  and  (i+2,j,k),  to  define  the  elements  of 
the  split  flux  Jacobian  matrices,  and  (2)  using  interpolation  with 
three  upwind  biased  points,  (i-l,j,k),  (i,j,k)  and  (i+l,j,k),  to 
define  Ui.  and  similarly  using  points  (i,j,k),  (i+l,j,k)  and 
(i+2,j,k)  to  define  Ur.  The  third  order  approximations  are 
important  within  boundary  layers  to  prevent  the  numerical 
dissipation  introduced  by  the  flux  splitting  from  significantly 
competing  with  viscous  terms  of  the  governing  differential 
e^ations.  Through  reflective  boundary  conditions,  third  order 
differencing  can  be  carried  all  the  way  to  the  wall.  Finally,  at 
shock  waves,  where  pressure  gradients  are  large,  and  the  flow  is 
supersonic  on  one  side  of  surface  i+l/2  and  subsonic  on  the  other, 
first  order  flux  splitting  is  used. 

For  each  mesh  point  in  a  three-dimensional  flow  the  algorithm 
uses  data  from  24  four  neighboring  points  to  evaluate  all  the 
special  derivatives  of  the  full  Navier-Stokes  equations  explicitly. 
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This  is  used  to  construct  the  explicit  driving  term,  sometimes 
called  the  residual  terra,  on  the  right  hand  side  of  the  difference 
equations.  The  left  hand  implicit  side  of  the  difference  equations 
uses  data  at  only  six  mesh  points  plus  the  central  point  itself. 
Only  the  thin  layer  Navier-Stokes  terms  are  evaluated  implicitly. 
The  mesh  points  are  located  at  the  centroids  of  finite  volumes 
forming  the  discretized  space  surrounding  the  body.  The  resulting 
block  seven-diagonal  matrix  equation  is  formulated  in  delta  law 
form.  That  is,  the  solution  change  per  time  step,  5U,  is  the 
unknown  to  be  solved  for.  The  matrix  equation  is  given  below. 


wheie 

1,J,K  ^2: 


The  hepta-diagonal  matrix  equation  is  inverted  by  relaxation 
as  follows.  The  mesh  points  are  ordered  into  an  i,j,k  array.  The 
values  for  the  unknown  iUs  are  initialized  by  setting  them  to  zero. 
For  each  i  plane  a  block  tridiagonal  inversion  in  the  j  direction 
is  carried  out  for  each  value  of  k  in  the  plane.  During  this 
inversion,  evaluations  needed  at  points  (i+l,j,k)  and  (i-l,j,k) 
adjacent  to  a  line  of  constant  i,k  use  the  latest  available  data  in 
line  Gauss-Seidel  fashion.  Evaluations  needed  at  off  line  points 
(i,j,k+l)  and  (i,j,k-l)  may  use  data  from  the  last  block 
tridiagonal  inversion  along  the  line  of  constant  i, j,  if  it  exists, 
in  line  Jacobi  fashion.  Yet,  in  the  present  code  an  axisymmetric 
approximation  for  these  terms  is  used  instead.  This  approximation 
effectively  replaces  these  two  terms  at  (i,j,k+l)  and  (i,j,k-l) 
with  a  new  one  at  the  diagonal  element  (i,j,k).  This  first  step 
procedure  is  shown  below  (see  top  of  the  next  page).  In  the 
equation  for  this  procedure  the  fiUs  at  (i,j+l.k),  (i,j,k)  and  (i,j- 
l,k)  are  solved  for  simultaneously. 

After  completion  of  the  first  step,  block  tridiagonal 
inversion  in  the  k  direction  proceeds  for  each  j  in  the  i  plane. 
As  before,  evaluations  at  (i+l,j,k)  and  (i-l,j,k)  are  made  in  line 
Gauss-Seidel  fashion  and  those  at  (i,j+l,k)  and  (i,j-l,k)  are  also 
made  in  Gauss-Seidel  fashion.  It  was  found  that  it  was  best  to 
order  these  k-direction  block  tridiagonal  Inversions  from  the  top 
down,  that  is,  in  the  descending  j  direction.  Thus,  the  newest 
information  is  moving  in  the  same  direction  as  the  principal 
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GAUSS-SEIDEL  LINE  RELAXATION  -  STEP  1 
BLOCK  TRIDIAGONAL  INVERSION  -  j  DIR. 


AS  ut}  ,^*Ss  + (?6  y;,"]. , ,  * + 

uf:).,,-,  ^GSuf;l,.,=  Ay,,,.,* 


{*)  indicates  latest  available  data 
(*)  indicates  axi symmetric  approximation 

component  of  the  flow  in  the  i  plane*  from  the  freestream,  through 
the  bow  shock  wave  and  towaid  the  body.  Therefore,  the  data  for  fiU 
at  (i,j+l,k)  comes  from  the  last  k-direction  line  inversion  at  j+1 
and  that  for  (i,  j“l,k)  comes  from  a  j-direction  line  inversion  from 
step  1.  This  second  step  is  given  below. 


GAUSS-SEIDEL  LINE  RETJiXATION  -  STEP  2 
BLOCK  TRIDIAGONAL  INVERSION  -  k  DIR. 


AS  y,',7,=-i  *§S  ^ 

ssu}:i,j_i,*ssuj‘ij,i,* 

utI%^*GS  yi,7,^i.,  =  Ay,,,,* 


(*)  indicates  latest  available  data 


After  the  inversion  procedures  for  both  the  j  and  k  directions 
within  plane  i  have  been  completed  the  next  i  plane  is  processed. 
In  general,  each  i  plane  is  processed  from  the  last  to  the  first 
and  then  a  second  time  from  the  first  back  to  the  last.  However, 
if  the  flow  is  supersonic  in  the  i-direction,  or  more  p'-ecisely,  if 
the  flow  component  normal  to  each  i  plane  is  supersonic  at  each 
mesh  point  in  the  flow,  except  within  the  boundary  layer  region, 
then  only  one  sweep  needs  be  made  from  the  inlet  i-plane  to  the 
exit  i-plane.  The  algorithm  requires  one  block  tridiagonal 
inversion  for  each  i,j  and  i,k  line  in  the  mesh  per  sweep  per  time 
step.  It  is  symmetric  with  respect  to  the  j  and  k  directions  but 
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not  the  i  direction. 


DISSIPATION 

There  is  no  added  dissipation  in  the  present  method.  But  this 
does  not  mean  that  there  is  no  dissipation  in  the  numerical  method. 
All  methods  require  numerical  dissipation  to  control  numerical 
Instabilities.  It  is  an  essential  ingredient.  The  present  method 
contains  dissipation  through  the  upwind  flux  split  difference 
approximations  to  the  Euler  terms  plus  the  centered  difference 
approximations  used  for  the  viscous  terms  only.  Dissipation 
introduced  this  way  is  said  to  be  genuine.  When  dissipation  is 
introduced  by  adding  new  terms  to  the  difference  equations  that  do 
not  correspond  to  terms  in  the  governing  differential  equations,  it 
is  artificial.  Some  justify  the  latter  kind  because  they  are  able 
to  control  the  amount  added  by  turning  up  or  down  a  parameter 
multiplying  the  added  term.  It  can  be  turned  up  to  the  point  of 
rendering  a  hyper-active  numerical  solution  comatose,  perhaps 
masking  the  effect  of  a  fundamental  bug  in  the  program.  If 
something  is  wrong  in  the  program,  it  is  better  to  have  a  quick 
death  than  a  lingering  program  life.  In  general,  humans  can  not  be 
trusted  with  artificial  viscosity  in  numerical  calculations. 

Yet  dissipation  is  essential.  The  successes  of  computational 
fluid  dynamics  can  be  explained  in  terms  of  the  numerical  methods 
having  the  right  dissipation.  Many  tried  unsuccessfully  to  solve 
the  transonic  small  disturbance  equation  using  central  differences 
at  each  mesh  point  before  Murman  and  Cole  [13]  made  the 
modification  to  use  backward  differences  in  supersonic  flow 
regions.  Central  differences  at  supersonic  flow  mesh  points  are 
anti-dissipative  and  backward  differences  are  dissipative.  Many 
rushed  to  use  this  strategy  on  the  full  potential  equation  without 
success  until  Jameson  [14]  devised  a  rotated  difference  scheme  that 
again  had  the  effect  of  replacing  anti-dissipative  difference 
approximations  with  dissipative  ones.  Neither  case  had  terms  added 
artificially  to  the  difference  equations. 

The  Steger-Warming  flux  vector  split  procedure  is  fairly 
dissipative,  too  dissipative  to  be  used  where  viscous  effects  are 
important.  It  was  originally  Intended  to  be  used  to  solve  the 
Inviscid  Euler  equations,  but  it  can  be  modified  to  greatly  reduce 
its  dissipation  so  that  it  can  be  used  for  viscous  flow  as  well. 
The  Roe  flux  difference  vector  splitting  procedure  can  also  be  used 
for  viscous  flow  as  it  is,  but  under  certain  flow  conditions 
additional  dissipation,  sometimes  called  entropy  correction,  needs 
to  be  added.  The  method  of  the  present  paper  uses  a  modified 
yersion  of  the  Steger-Warming  procedure  described  in  Ref. [19].  It 
is  used  throughout  the  flow  field  in  both  viscous  and  inviscid  flow 
regions.  However,  it  is  not  always  sufficient  to  control  numerical 
difficulties  and  in  some  regions  of  the  flow  the  original  more 
dissipative  Steger-Warming  procedure  is  blended  in  with  it.  These 
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difficulties  occur  at  shock  waves. 


DIFFICULTIES  AT  SHOCKS 

The  numerical  difficulties  at  shock  waves  can  be  categorized 
as  either  annoying  or  catastrophic.  The  annoying  consist  of 
numerical  shock  buzz,  which  is  the  continual  moving  of  the  shock 
forward  and  backwards  without  ever  settling  down.  This  motion  can 
be  either  periodic  or  a  chaotic  non-repeating  motion.  When  it  is 
present  it  occurs  near  the  stagnation  streamline  of  the  flow.  It 
feeds  on  a  supersonic-subsonic  interaction  at  the  bow  shock  wave. 
The  catastrophic  consist  of  negative  temperatures  and  densities 
occurring  at  the  foot  of  the  shock,  also  usually  located  near  the 
stagnation  streamline.  These  difficulties  are  magnified  as  the  CFL 
number  is  increased.  For  this  study  CFL  numbers  of  the  order  of 
one  million  were  used.  The  workstation  word  size  essentially 
limited  the  CFL  to  this  value,  not  the  numerical  method. 

The  alignment  of  the  grid  with  the  shock  near  the  stagnation 
streamline  is  a  key  factor  in  producing  or  eliminating  these 
difficulties.  Perfect  alignment  along  the  entire  shock  is 
impractical  unless  the  mesh  is  actively  shock-fit  during  the 
calculation.  It  was  observed  that  where  errors  in  the  alignment 
occur,  it  was  better  for  the  shock  to  cut  through  the  mesh  toward 
the  body  rather  than  toward  the  freestream.  This  tends  to  reduce 
the  feedback  or  interaction  mechanism  within  the  subsonic  region 
behind  the  strong  part  of  the  bow  shock. 

Secondly,  as  mentioned  above,  pure  Steger-Warming  was  blended 
in  with  modified  Steger-Warming  flux  splitting  at  shock  wave 
locations.  The  amount  used  depended  upon  the  pressure  gradient, 
PG.  The  weight  function  used  for  it  varied  from  zero  to  one 
according  to  the  formula 

WT  =  1.0-1.0/(1.0+PG*PG) 

Analysis  of  the  Steger-Warming  algorithm  shows  that  there  is  a 
considerable  mixing  of  fluid  across  mesh  surfaces.  This  should 
lead  to  a  smearing  of  the  shock  that  prevents  negative  densities 
and  temperatures.  This  is  exactly  what  occurs  when  it  is  used 
explicitly.  Densities  and  energies,  both  positive  quantities,  are 
mixed.  However,  when  it  is  also  used  implicitly  as  in  the  delta 
law  form  of  the  difference  equations  given  above,  it  mixes  the 
"changes'*  in  densities  and  energies.  These  changes  can  be  of 
either  sign.  Negative  changes  can  be  mixed  forward  through  the 
shock  wave  causing  negative  densities  and  energies  in  the 
freestream  at  the  foot  of  the  shock.  To  avoid  this  from  happening, 
the  blending  in  of  pure  Steger-Warming  was  used  only  to  obtain  the 
explicit  in  the  above  equations  and  was  not  used  at  all  on 
the  Implicit  side,  the  left  side,  of  the  block  matrix  equations 
given  above. 
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STABILITY 


The  present  numerical  technique  is  implicit  and  can  be  run 
with  CFL  numbers  as  high  as  a  billion.  Yet,  a  calculation  cannot 
be  started  with  any  arbitrarily  large  time  step  size.  Initially 
the  transients  in  the  solution  are  large  and,  even  though  the 
method  may  be  stable  according  to  linear  theory  for  any  time  step 
size,  it  is  limited  in  practice  because  of  the  non-linearity  of  the 
governing  equations  being  solved.  One  would  expect  this  if  the 
Jacobian  matrices  appearing  on  the  implicit  side  of  the 
approximating  finite  difference  equations,  or  as  in  the  present 
case  finite  volume  equations,  change  significantly  during  a  single 
time  step.  In  impulsively  started  viscous  flow  calculations, 
exceptionally  high  shear  exists  at  no  slip  boundary  surfaces. 
Large  transients  can  also  occur  in  the  flow  field  after  the  start 
of  the  calculation,  perhaps  associated  with  the  movement  of  a 
strong  bow  shock  wave.  In  a  large  three-dimensional  calculation 
it  is  difficult  to  anticipate  all  such  transients  and  then  to 
determine  how  large  the  time  step  can  be  a  priori.  Each 
calculation  can  represent  a  significant  investment  in  computer  time 
which  can  be  lost  if  it  becomes  unstable. 

To  avoid  this  difficulty,  a  new  subroutine  was  built  into  the 
present  three-dimensional  Navier-Stokes  program  to  automatically 
adjust  the  time  step  through  periods  of  strong  flow  transients.  A 
time  step  is  chosen  initially  rather  arbitrarily  and  increased 
geometrically  during  the  calculation.  At  each  time  step,  after  the 
changes  in  the  solution  have  been  computed  implicitly,  the  changes 
are  checked  to  determine  if  either  density  at  any  mesh  point  in  the 
flow  will  change  by  more  than  a  factor  of  one-half,  or  if  an 
internal  energy  any  where  in  the  flow  field  will  decrease  by  more 
than  a  factor  of  one-half.  If  either  condition  is  true,  the  time 
step  during  that  step  is  cut  to  that  value  for  which  either  the 
density  or  internal  energy  will  be  limited  to  a  change  of  at  most 
half.  The  solution  changes  do  not  have  to  be  recalculated 
implicitly  because  they  depend  linearly  on  the  time  step  size. 
They  need  only  to  be  cut  everywhere  by  the  same  factor  the  time 
step  is  cut.  An  alternative  version  of  cutting  the  time  step 
locally  was  tried  with  some  success  also,  but  it  introduces 
non-conservation  into  the  calculation  which  was  not  preferred.  For 
the  computed  results  to  be  presented,  reductions  in  time  step  size 
were  observed  periodically  during  the  early  part  of  the 
calculations,  but  then  the  need  for  time  step  reductions 
disappeared  entirely  as  convergence  occurred.  In  addition,  the 
time  step  was  increased  continuously  during  the  calculation  at  a 
rate  that  doubled  it  every  eight  steps. 


COMPUTATIONAL  RESULTS 

The  computer  code  was  applied  to  two  test  flow  problems  for 
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which  experimental  or  flight  data  was  available  for  program 
validation. 

Lobb  Test  Case 

The  first  case  simulated  the  experiment  of  Lobb  in  1964  [1]  for 
Mach  15.3  flow  past  a  sphere  of  radius  0.635cm.  and  Reynolds  number 
1.47x10*  in  air.  The  freestream  pressure  was  664  Netrtons  per 
square  meter,  the  freestream  temperature  was  293K  and  the  body 
surface  temperature  was  fixed  at  lOOOK.  The  shock  wave  position 
was  guessed  at  the  start  of  the  calculation  and  inviscid  shock  wave 
theory  was  used  to  define  the  initial  solution.  Two  cases  were 
run,  perfect  gas  and  seven  species  reacting  air,  on  an  axisymmetric 
24x30  mesh  shown  in  Fig.l.  Care  was  taken  to  fit  the  mesh  to  the 
shock  near  the  axis  of  symmetry.  The  shock  wave  locations  for  each 
case  are  shown  in  Fig. 2  compared  with  the  experimental  data  of  Lobb 
obtained  from  a  Schlieren  photograph.  The  shock  wave  represents  in 
a  sense  the  integrated  effect  of  the  state  of  the  gas  between  the 
body  and  shock.  Note  the  nonequilibrium  results  compare  much 
better  with  the  experimental  results,  indicating  significant 
nonequilibrium  flow  is  present.  Both  Candler  [3]  and  Gokcen  [4] 
used  this  test  case  for  code  validation  with  excellent  agreement 
with  experiment.  Notice  the  shock  ripples  near  the  exit  as  the 
shock  crosses  through  the  mesh.  Although  the  shock  simulation  is 
degraded  in  this  region,  it  did  not  cause  any  of  the  difficulties 
discussed  earlier  because  it  occurred  dotmstream  of  the  sonic  line. 

The  pressure  contours  for  nonequilibrium  flow  are  shown  in 
Fig. 3  and  mass  fractions  for  specie  NO  are  shown  in  Fig. 4  Note 
that  the  NO  mass  fractions  concentrate  in  a  region  between  the  body 
and  shock.  The  mass  fractions  of  NO  are  zero  in  the  freestream, 
build  up  upon  passing  through  the  shock  wave  where  dissociated 
nitrogen  and  oxygen  are  found  and  then  decrease  as  NO  dissociates 
in  the  wall  cooled  boundary  layer.  The  wall  was  assumed  to  be 
noncatalytic.  The  stagnation  streamline  temperatures  are  shown  in 
Fig. 5.  In  the  figure  the  body  is  located  at  x=0  and  the  shock 
approximately  x/r=-0.085.  Note  that  the  rotational  temperature  is 
in  equilibrium  with  the  translational  temperature  everywhere  except 
near  the  shock,  and  the  vibrational  temperature  is  out  of 
equilibrium  everywhere  except  where  the  stagnation  streamline 
approaches  the  body  fixed  wall  temperature  of  lOOOK.  Figures  6-8 
show  contours  of  translational,  rotational,  and  vibrational 
temperatures.  Note  again  that  the  rotational  temperature  is  in 
equilibrium  with  the  translational  temperature  almost  everywhere 
and  that  the  vibrational  temperature  is  out  of  equilibrium  almost 
everywhere.  These  results  agree  with  those  of  Candler  [3]  who 
assumed  rotational  equilibrium  and  with  Gokcen  [4]  who  did  not. 

RAM-C  II  Flight  Measurements 

During  the  late  1960 's  sphere  cone  bodies  were  rocket 
launched,  turned  around  at  high  altitudes  and  fired  down  through 
the  upper  atmosphere  [ 2 ] .  Microwave  antenna  carried  on  board  the 
vehicle  made  measurements  of  peak  electron  number  densities.  This 
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flow  field  was  simulated  numerically  by  Candler  [3]  with  excellent 
agreement  at  all  altitudes  calculated.  In  the  present  study,  only 
one  altitude  at  71  km  will  be  considered.  The  sphere-cone  body  of 
nose  radius  0.1524m,  cone  half  angle  of  9°  and  body  length  of 
1.295m  traveled  at  Mach  25.9  The  freestreara  pressure  was  4.898 
Ne%rtons  per  square  meter  and  the  freestream  tei*>perature  was  216K. 
The  body  surface  was  assumed  to  be  noncatalytic  and  at  a  fixed 
1500K.  The  axisymmetric  24x30  mesh  and  pressure  contours  are  shovm 
in  Figs.  9  and  10. 

The  translational,  rotational  and  vibrational  temperature 
profiles  along  the  stagnation  streamline  are  shown  in  Fig. 11.  This 
flow  exhibited  considerable  translational,  rotational  and 
vibrational  nonequilibrium  before  arriving  at  the  fixed  wall 
temperature.  The  shock  wave  location  appears  to  be  approximately 
x/r=-0.1.  Figure  12  shows  the  mass  fractions  of  Nj,  Oj,  NO,  N  and 
O  also  along  the  stagnation  streamline.  Note  that  the  shock  wave 
appears  here  to  be  located  at  only  x/r=-0.07.  The  difference  is 
explained  by  the  need  to  build  up  a  significant  vibrational 
temperature  before  dissociation  of  N^  and  Oj  can  occur  as  is 
formulated  by  Park's  geometric  mean  temperature  [3],  (T*T„ib)^^*» 
used  in  the  chemical  rate  equations.  Note  also  that  the  species 
nearly  all  return  to  their  freestream  values  at  the  wall.  Although 
the  wall  is  noncatalytic,  it  is  sufficiently  cool  to  bring  about 
recombination  of  N,  and  O,.  The  contours  of  N,  and  O,  are  given  in 
Figs.  13  and  14.  Note  also  here  that  the  species  return  to  near 
their  freestream  values  along  the  whole  length  of  the  body. 
Figures  15-17  show  contours  of  translational,  rotational  and 
vibrational  temperatures.  Note  that  although  these  three 
temperatures  were  far  out  of  equilibrium  along  the  stagnation 
streamline,  the  rotational  temperature  is  in  fairly  good  agreement, 
as  assumed  by  Candler  [3],  for  most  of  the  flow  field  and 
vibrational  temperature  remains  in  nonequilibrium  almost 
everywhere . 

Figure  18  compares  the  measured  peak  electron  number  densities 
with  the  computation.  The  flow  field  was  calculated  in  two  parts, 
the  nose  srvrtion,  x  <  2.4,  and  the  section  downstream  of  the  nose, 
X  >  2.4.  Ajl though  this  flow  could  have  been  calculated  on  a  single 
mesh,  this  test  served  to  illustrate  the  sectional  calculation  of 
a  supersonic  flow  field.  This  is  more  efficient  because  no 
calculations  need  be  made  until  the  flow  section  immediately  ahead 
is  converged.  This  figure  also  illustrates  the  speed  of 
convergence  of  the  present  numerical  method.  After  the  nose 
section  converged,  the  data  at  the  nose  section  exit  was  used  to 
initialize  the  flow  field  of  the  downstream  section.  Convergence 
was  then  obtained  in  only  about  16  time  steps.  Note  the  comparison 
with  the  experiment  is  not  perfect,  but  also  note  that  there  is  a 
four  decade  variation  in  number  densities  sho«m  within  this  figure. 
The  thought  of  getting  as  close  to  the  measured  values  as  shorn 
here  was  only  a  hope  a  priori  to  performing  the  calculation.  For 
this  same  altitude  of  71  km.,  Candler  [3]  fit  the  data 
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exceptionally  well  on  a  finer  grid. 

Generic  Hypersonic  Vehicle 

The  computer  program  was  applied  to  simulate  the  flow  about  a 
generic  hypersonic  vehicle  at  Mach  25  and  at  high  altitude.  The 
body  consisted  of  a  sphere-cone-cylinder  with  a  delta  wing 
attached.  The  delta  wing  tips  were  folded  up  to  resemble  a 
hypersonic  glider.  An  engine  was  then  placed  below  and  attached  to 
the  wing-fuselage.  The  underside  of  the  body  downstream  of  the 
engine  exit  was  contoured  to  form  a  nozzle  surface.  A  side  view  of 
the  body  and  its  dimensions  are  given  in  Fig.  19  and  flight 
conditions  are  given  in  Fig. 20  with  a  rotated  side  view  of  the 
vehicle.  The  computational  strategy  is  given  in  Fig. 21  with  a  view 
of  the  underneath  side  of  the  vehicle.  The  flow  was  cut  into  four 
sections,  nose,  fore  body,  engine  and  after  body  nozzle  sections. 
Because  of  body  and  flow  symmetry,  only  the  flow  on  one  side  of  the 
body  symmetry  plane  needed  to  be  computed.  Each  section  was  run  in 
order  from  the  nose  to  the  tail,  and  only  after  the  preceding 
section  converged.  The  solution  at  the  exit  plane  of  each  section 
was  used  to  initialize  the  flow  in  the  section  just  downstream  of 
itself.  If  the  geometry  at  the  exit  resembled  the  geometry  of  the 
neighboring  dotmstream  section,  as  is  the  case  near  the  nose,  this 
initialization  facilitated  fairly  rapid  numerical  convergence, 
otherwise,  as  at  the  last  section  of  the  body,  convergence  was 
slower. 

Heat  transfer  contours  on  the  nose  surface  are  shown  in 
Fig. 22.  Care  must  be  used  in  reading  values  from  the  key  or  table 
at  the  right  side  of  the  figure  because  the  highest  values  shown  in 
kilowatts  per  square  meter  have  been  divided  by  10  by  the  plotting 
software.  Stagnation  point  heating  is  4.9x10'^  and  heating  along 
the  leading  edge  of  the  emerging  wing  approaches  2.8x10''  kilowatts 
per  square  meter  at  the  exit  of  the  nose  section. 

A  surface  grid  for  the  fore  body  section  is  shown  in  Fig. 23 
and  a  cross  section  of  the  mesh  along  the  symmetry  plane  is  shown 
in  Fig. 24.  The  mesh  has  thirty  finite  volume  cells  from  body  to 
free  stream.  The  flow  entering  the  engine  inlet  plane  was 
accelerated  to  produce  thrust  at  the  exit  plane  without  changing 
the  mass  flow  rate.  The  temperatures  of  the  gas  entering  the 
engine  inlet  was  approximately  the  same  as  that  leaving  the  engine 
exit.  The  flow  acceleration  was  designed  to  be  approximately  equal 
to  the  sum  of  the  vehicle  surface  drag  and  engine  ram  drag.  The 
flow  at  each  mesh  point  leaving  the  fore  body  mesh  and  entering  the 
engine  inlet  was  accelerated  and  mapped  one  to  one  with  mesh  points 
entering  the  after  body  nozzle  section  from  the  nozzle  exit.  Other 
values  at  the  inlet  planes  of  the  engine  and  after  body  nozzle 
sections  were  obtained  by  interpolation  from  the  corresponding  exit 
planes.  Pressure  contours  and  velocity  vectors  in  the  syrwetry 
plane  are  shown  in  Figs.  25  and  26.  Note  the  clearly  defined 
laminar  boundary  layers  appearing  ir  ’*.he  velocity  vector  plot. 
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The  surface  mesh  of  the  generic  hypersonic  vehicle  is  shown  in 
Fig. 27.  The  trailing  edge  of  the  vehicle  is  a  knife  edge  as  are 
the  engine  inlet  and  exit  surfaces.  Surface  pressure  contours  in 
Newtons  per  square  meter  are  shorn  in  Fig. 28.  The  peak  surface 
pressure  occurs  at  the  stagnation  point  and  achieves  a  value  of 
4.1x10^.  High  surface  pressures  also  occur  along  the  leading  edges 
of  the  folded  delta  wing.  Surface  heat  transfer  contours  are  sho«fn 
in  Fig. 29  again  with  peak  heating  at  the  nose  and  wing  leading 
edges.  Total  thrust,  lift,  drag  and  heat  transfer  for  the  generic 
hypersonic  vehicle  are  given  in  Fig. 30.  Engine  thrust  was  adjusted 
to  give  a  net  forward  impulse  to  the  vehicle. 

A  summary  of  the  bookkeeping  in  terms  of  section  mesh  sizes 
and  number  of  time  steps  to  convergence  is  given  in  Fig. 31. 
Section  4,  the  after  body  nozzle  section,  took  the  longest  to 
converge.  The  flow  in  this  section  had  the  longest  way  to  go  from 
the  initial  condition  to  converged  solution,  because  of  extreme 
flow  conditions  in  the  expanding  nozzle  region  downstream  of  the 
engine  exit.  The  automatic  time  step  reduction  procedure, 
described  in  the  stability  section  of  this  paper,  was  called  at 
each  time  step  during  the  first  32  time  step  run,  thus  reducing  the 
actual  flow  time  simulated.  The  program  was  run  an  additional  64 
time  steps,  the  last  half  of  which  did  not  require  automatic  time 
step  reductions.  At  approximately  96  time  steps,  the  flow  settled 
out. 


The  computational  cost  of  this  three-dimensional  Gauss-Seidel 
line  relaxation  method  on  a  mesh  of  size  ILxJLxKL  is  approximately 
the  time  taken  by  ILx(JL+KL)  block  tridiagonal  inversions.  On  a 
SUN  SPARCstation  1,  a  section  mesh  of  24x31x26,  for  32  time  steps, 
integrating  the  13  equations  governing  chemically  reacting  thermal 
nonequilibrium  flow,  required  approximately  5  days  to  finish.  The 
flow  field  calculation  for  each  section  was  run  in  the  background 
while  the  workstation  was  used  also  to  set  up  the  mesh  and  initial 
conditions  required  for  the  next  flow  section.  However,  the  human 
times  required  to  set  up  the  surface  geometry  and  mesh  for  the  next 
section  to  be  done  were  usually  longer  than  5  days.  Thus,  the 
SPARCstation  1,  relatively  slow  by  today's  standards,  was  sitting 
idle  much  of  the  time.  Nevertheless,  this  calculation  of  the 
'nonequilibrium  flow  past  a  complete  hypersonic  vehicle  in  powered 
flight  is  believed  to  be  the  first  anywhere,  even  at  aeronautical 
laboratories  with  supercomputer  access. 


CONCLUSION 

An  efficient  numerical  method  has  been  presented  for  solving 
the  compressible  Navier-Stokes  equations  in  three  dimensions.  It 
has  been  extended  to  solve  for  the  chemically  reacting  flow  in 
thermal  nonequilibrium  encountered  under  hypersonic  conditions. 
Experimental  and  flight  test  data  have  been  used  to  validate  a 
computer  program  based  upon  this  method .  The  program  has  been 
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used,  on  an  inexpensive  computer  workstation,  to  predict  the 
complete  flow  field  about  a  generic  hypersonic  vehicle  in  powered 
flight. 
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Fig.l  Computational  mesh  for  Lobb  test  case. 


Fig. 2  Comparison  of  shock  wave  locations,  Lobb  test  case. 
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Fig. 3  Pressure  contours,  Lobb  test  case. 


Fig. 4  NO  mass  fractions,  Lobb  test  case. 
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Fig. 5  Stagnation  streamline  temperature  profiles,  Lobb  test  case. 


Print  II  lobbne46t.plt  ||  lobbne4a.kpln.t«mp 

TRANSLATIONAL  TEMPERATURE  CONTOURS 


Fig. 6  Translational  temperature  contours,  Lobb  test  case. 
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Fig. 7  Rotational  temperature  contours,  Lobb  test  case. 


print  II  lobbnp48t.p(t  H  lobbn^6.kptn.tqmp 


VIBRATIONAL  TEMPERATURE  CONTOURS 


Fig. 8  Vibrational  temperature  contours,  Lobb  test  case. 
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Fig. 9  Computational  Mesh,  RAM~C  II  test  case. 


Fig. 10  Pressure  contours,  RAH-C  II  test  case. 
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Fig. 11  Stagnation  streamline  temperature  profiles,  RAM-C  II  test 
case. 


Fig. 12  Stagnation  streamline  mass  fractions,  RAM-C  II  test  case. 
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PEAK  ELECTRON  NUMBER  DENSITY 

RAM-C  TEST  CASE.  71km.  altitude 
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Fig.  18  Peak  electron  nustber  density,  RAM-C  II  test  case. 
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APPLICATION  TO  A  GENERIC  HYPERSONIC  VEHICLE 


SPHERE-CONE-CYLINDER-DELTA  WING-BODY 
LENGTH  =  71m..  SPHERE  RADIUS  =  1  m. 
CYLINDER  RADIUS  =  4m. 

CONE  HALF  ANGLE  =  7.5  DEGREES 
DELTA  WING  HALF  ANGLE  =  20  DEGREES 


Fig. 19  Generic  hypersonic  vehicle  definition. 


pO)  Print  II  ghvasff.pit  ||  srfdat.pdat 


FLIGHT  CONDITIONS 


LEVEL  FLIGHT  AT  71km.  ALTITUDE 
MACH  25  LAMINAR  FLOW 
FREESTREAM  PRESSURE  4.9  N/m**2 
FREESTREAM  TEMPERATURE  216  DEGREES  K 
BODY  TEMPERATURE  1500  DEGREES  K 


Fig. 20  Flight  conditions  for  generic  hypersonic  vehicle. 
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Fig. 21  Computational  strategy  for  hypersonic  vehicle  flow  field. 


Fig. 22  Heat  transfer  at  nose  section  of  hypersonic  vehicle. 
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Fig. 23  Surface  nesh  of  fore  body  section  of  hypersonic  vehicle. 


Print  li  ghvas.pit  H  ghvaside _ 

,  SYMMETRY  PLANE  MESH 

GENERIC  HYPERSONIC  VEHICLE 
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Fig. 25  Pressure  contours  at  syaaetry  plane  of  hypersonic  vehicle. 
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Fig, 27  Surface  mesh  of  hypersonic  vehicle. 


Fig. 28  Surface  pressure  contours  of  hypersonic  vehicle 


[30)  Print  II  ntivaal.plt  ||  trtdat.p55~ 
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Fig. 29  Heat  transfer  contours  at  surface  of  hypersonic  vehicle. 


II  ghva-plt  II  Qfwl.aide 

THRUST.  DRAG  AND  HEAT  TRANSFER 

GENERIC  HYPERSONIC  VEHICLE 

SURFACE  AREA  -  4,136  square  meters 
TOTAL  HEAT  TRANSFER  -  68,262  kwatts 

LIFT  -  10,333  Newtons 

SURFACE  DRAG  =  46,741  Newtons 

ENGINE  RAM  DRAG  =  189,328  Newtons 

TOTAL  DRAG  =  236,069  Newtons 

ENGINE  THRUST  =  250,140  Newtons 

Fig. 30  Thrust,  lift,  drag  and  heat  transfer  for  hypersonic  vehicle. 
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Print  II  ghvaplt  ||gh»1.tide 

BOOKKEEPING 

GENERIC  HYPERSONIC  VEHICLE 

SECTION  1  -  NOSE  SECTION 

MESH  -  1 6X31 X26.  32  TIME  STEPS 

SECTION  2  -  FORE  BODY  SECTION 
MESH  -  24X31 X26.  32  TIME  STEPS 

SECTION  3  -  ENGINE  SECTION 

MESH  -  24X31X26,  40  TIME  STEPS 

SECTION  4  -  AFTER  BODY/NOZZLE  SECTION 
MESH  -  24X31X26,  96  TIME  STEPS 


Fig. 31(a)  Section  mesh  sizes  and  time  steps  to  convergence  for 
numerical  solution  of  flow  past  a  generic  hypersonic 
vehicle  in  powered  flight. 


;3D)  Print  II  ghvasrt.pit  ||  ghvasri.pdat 

BOOKKEEPING  CONTINUED 

GENERIC  HYPERSONIC  VEHICLE 


COST  OF  3-D  GAUSS-SEIDEL  LINE  RELAXATION 

ON  A  MESH  OF  DIMENSION  IL  x  JL  x  KL 
IL(JL+KL)  BLOCK  TRIDIAGONAL  INVERSIONS 
ARE  REQUIRED  PER  TIME  STEP 

ROUGHLY  ONE  WEEK  OF  COMPUTING  TIME 
PER  SECTION  ON  A  SUN  SPARCstation  1 


Fig. 31(b)  Computational  cost  of  numerical  solution  for  flow  past 
generic  hypersonic  vehicle  in  powered  flight. 
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